Towards adiabatic waveforms for inspiral into Kerr black holes: 
II. Dynamical sources and generic orbits 
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This is the second in a series of papers whose aim is to generate "adiabatic" gravitational wave- 
forms from the inspiral of stellar-mass compact objects into massive black holes. In earlier work, we 
presented an accurate (2+l)D finite-difference time-domain code to solve the Teukolsky equation, 
which evolves curvature perturbations near rotating (Kerr) black holes. The key new ingredient 
there was a simple but accurate model of the singular source term based on a discrete represen- 
tation of the Dirac-delta function and its derivatives. Our earlier work was intended as a proof of 
concept, using simple circular, equatorial geodesic orbits as a testbed. Such a source is effectively 
static, in that the smaller body remains at the same coordinate radius and orbital inclination over 
an orbit. (It of course moves through axial angle, but we separate that degree of freedom from the 
problem. Our numerical grid has only radial, polar, and time coordinates.) We now extend the 
time-domain code so that it can accommodate dynamic sources that move on a variety of physi- 
cally interesting world lines. We validate the code with extensive comparison to frequency-domain 
waveforms for cases in which the source moves along generic (inclined and eccentric) bound geodesic 
orbits. We also demonstrate the ability of the time-domain code to accommodate sources moving on 
interesting non-geodesic worldlines. We do this by computing the waveform produced by a test mass 
following a "kludged" inspiral trajectory, made of bound geodesic segments driven toward merger 
by an approximate radiation loss formula. 

PACS numbers: 04.25.Nx, 04.30.Db, 04.30.-w 
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I. INTRODUCTION 



A. Background 



The extreme mass ratio limit of general relativity's 
two-body problem has been a major focus of work in 
recent years. This limit corresponds to a stellar mass 
compact object that orbits and perturbs a massive black 
hole. The system generates gravitational waves (GWs) 
which drive the small body to inspiral into the large 
black hole. Measuring such "extreme mass ratio inspi- 
ral," or EMRI, events is a major goal for space-based 
GW antennae, particularly the LISA mission 1 . EMRIs 
should be measurable to a redshift z ~ 0.5 — 1. The 
event rate at this range is estimated to be high enough 
that a multiyear LISA mission should measure dozens 
to hundreds of EMRI events Because the smaller 
body only slightly perturbs the larger black hole's space- 
time, EMRI GWs are expected to provide an exception- 
ally clean probe of black hole properties. We expect to 
use EMRIs to measure black hole masses and spins with 
extremely good accuracy @], and even to test how well 
the spacetime meets the rather stringent constraints that 
the "no-hair" theorems of general relativity impose on 
black holes 



1 http://lisa.nasa.gov, http://sci.esa.int/lisa 



Understanding EMRI sources will require us to com- 
pare measured waves with theoretical models that are 
as accurate as possible. This goal motivates much re- 
cent EMRI work. The waves are sufficiently complicated 
that simply detecting them in LISA's datastream will be 
a challenge. Techniques for finding these events are cur- 
rently being developed and tested through the "Mock 
LISA Data Challenges" , or MLDCs (see Refs. @, 1 for 
overviews of recent MLDCs). An important input to 
these challenges (and to the development of EMRI mea- 
surement techniques more generally) are waveform mod- 
els that capture the true complexity of EMRI events (see 
[H, |T(| for discussion of recent work to include EMRI 
waves in the MLDCs). 

This paper presents a further step in our program to 
construct accurate EMRI wave models. As discussed in 
the introduction to Ref. [ill (Paper I), our goal is to 
make "adiabatic" waveforms — waveforms built by sep- 
arately treating the long-time dissipative evolution and 
the short-time conservative motion. In our present anal- 
ysis, we take the short-time motion to be a geodesic orbit 
of the background spacetime; our approach thus amounts 
to approximating the inspiral trajectory as a sequence of 
geodesic orbits. As discussed by Pound and Poisson 
this limit is more properly a "radiative" or "dissipative" 
approximation, since we do not include conservative self 
interactions. It may be possible to augment this analy- 
sis with at least some conservative effects [13| , so we be- 
lieve the program we are developing is capable of building 
truly adiabatic inspiral waveforms as described in 
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We will describe our goal as "adiabatic" waveforms, but 
the reader should bear in mind that the approximation 
we are currently developing is more restricted than this. 

Geodesic orbits are described (up to initial conditions) 
by three conserved constants: energy E, axial angular 
momentum L z , and "Carter constant" Q. Using black 
hole perturbation theory, we compute the rate at which 
these three constants evolve; fast and accurate frequency- 
domain codes make it possible to compute these rates 
of change fairly easily [3, HH EBj- We then build the 
parameter-space trajectory [E(t), L z (t),Q(t)] followed by 
the small body; choosing initial conditions, it is simple 
to build the coordinate-space worldline [r(t), 9(t), 4>(t)] of 
a particular inspiral. From this worldline, we build the 
source to a time-domain code. The output of this code 
is, at last, our model EMRI wave. 



B. Time-domain black hole perturbation theory 

Since the frequency-domain portion of this program 
is already well in hand, our current focus is on the time- 
domain code. In essence, our goal is to build a code which 
takes as input any physically reasonable worldline, and 
provides as output the waveform produced by a small 
body on this worldline. In Paper I, we demonstrated an 
accurate (2+l)D numerical code to solve, in the time do- 
main, the wave equation for curvature perturbations to 
a black hole — the Teukolsky equation 17]. Our code 
evolves the Weyl curvature scalar ^4, constructed by pro- 
jecting the vacuum curvature onto appropriate compo- 
nents of a null tetrad; see Paper I for details. The az- 
imuthal dependence of ^4 is separated out (due to the <p 
symmetry of black holes); the dependence on the Boyer- 
Lindquist coordinates r, 9, and t is found by evolving $4 
on a numerical r-9 grid. 

As is common in black hole perturbation theory, we 
treat the smaller body as a Dirac-delta point particle, 
leading to a singular source for the Teukolsky equation. 
In the frequency domain, the delta can be dealt with 
analytically, and presents no great challenge. By con- 
strast, accurately computing the effect of a sharp source 
on the time-domain code's numerical grid can be ex- 
tremely challenging. In Paper I, we presented a new tech- 
nique for treating the singular source term. Our innova- 
tion was to model the delta as a series of finite impulses, 
with the largest impulse located close to the delta's ar- 
gument, falling off rapidly as we move away from this 
"central" spike. Importantly, this approach allows us to 
accurately model the derivatives of the delta function. 
Since the Teukolsky equation source depends on first and 
second derivatives of the delta (as well as the delta it- 
self), this appears to give us an accuracy boost relative 
to other finite-difference delta representations (such as a 
truncated Gaussian), which may accurately capture the 
delta's behavior, but not do so well with the derivatives. 



C. This paper 

Paper I focused on the properties of this new source 
representation. To clarify this focus, we studied very 
simple orbits: We only considered the (astrophysically 
unlikely) case of circular, equatorial black hole orbits. 
We now extend this to include inclined, eccentric and 
generic orbits, as well as non-geodesic inspiral sequences. 

A particle in a circular, equatorial orbit has constant 
radial and angular coordinate, confining it to a fixed lo- 
cation on the r-9 grid. Eccentricity means that the orbit 
oscillates radially, crossing radial grid zones. Similarly, 
orbital inclination results in angular grid crossing. We 
quickly discovered that these new motions introduce high 
frequency numerical noise. This noise can be controlled 
by combining a low pass filter with a higher order dis- 
cretization of the delta function; details are given in Sec. 
ITT1 Aside from this mild extension of the basic formal- 
ism presented in Paper I, it was not terribly difficult to 
use our new source term to handle a broad class of as- 
trophysically interesting orbits. To validate our results, 
we present in Sec. Mil extensive comparisions with wave- 
form snapshots computed in the frequency domain [lil ]. 
demonstrating graphically and quantitatively (with ap- 
propriate overlap integrals) excellent agreement between 
the two techniques. 

As extensively discussed in the introduction to Paper 
I and here, our goal is to compute the waves from inspi- 
ral of a small body through a sequence of orbits. As a 
proof-of-concept demonstration of the feasibility of this 
idea, we present a simple example of inspiral in Sec. IIVI 
In this example, we evolve through our geodesic sequence 
using a "kludged" approximation to the rates of change 
of orbital constants, using the code described in Refs. 
[HI, [li| . These waveforms are not reliable EMRI models, 
but they illustrate the ease with which we can handle the 
effect of radiation emission on the motion of the source. 
Computing waves from an inspiral is no more of a com- 
putational challenge than computing waves from a bound 
geodesic. 

The next step will be to combine accurate radiative 
backreaction with our time-domain solver to compute 
"adiabatic" EMRI waveforms (albeit ones that still ne- 
glect conservative self interactions). Plans for this next 
step are described in our final summary, Sec. |Vj 



II. DYNAMICALLY VARYING DISCRETE 
DELTA FUNCTIONS 

In Paper I, we presented a method for representing 
a Dirac delta function and its derivatives on a discrete 
numerical grid. In that paper, we only considered a delta 
with fixed radial and angular position. Naive application 
of the discrete delta models presented in Paper I leads 
to instabilities when the particle moves in the numerical 
grid. The following argument outlines the root cause 
of these instabilities. Consider the function 8[x — a(t)], 
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where Xk < a(t) < Xk+i', i-e., the delta's peak varies 
with time and lies between two discrete grid points. Let 
Xi represent any discrete point on our grid, and let h = 
Xk+i — Xk = Xk — Xk-i be the grid resolution. Naive 
application of the results from Paper I might lead us to 
model the delta function with the impulse weights 

Si(tn) = a{tn \ ~ Xh fori = fc + l (2.1) 



h 2 

Xk+i - a(t n ) 
h 2 

everywhere else 



for i = k 



(2.2) 
(2.3) 



(This "two impulse" delta is in fact just the simplest 
representation we developed in Paper I, but is useful for 
the following discussion.) Each t n defines a time slice of 
our r-6 grid. As a varies from one time slice to another, 
so do the coefficients at Xk and Xk+i- The frequency 
spectrum of Sk(t n ) and dk+i(t n ) will reflect the amount 
of variation in a. A large variation in a will produce 
a high frequency component in the Fourier transform of 
the time series of each weight. These variations couple 
to the time derivatives in the homogeneous part of the 
Teukolsky equation. Consequently, the solution contains 
spurious high frequency features of numerical origin. 

Consider the extreme limit of this effect: a changes so 
rapidly that the delta's peak moves across a grid zone in 
a single time slice: 



so that 



but 



a(ti) = a(to) — h 



Xk < a(to) < Xk+i 



Xk- 



-l < a{ti) < x k 



(2.4) 



(2.5) 



(2.6) 



The weight of the delta function very suddenly becomes 
zero at Xk+i as we step from t = to to t = t\; likewise, 
the weight at Xk-i very suddenly jumps from non-zero to 
zero in this step. The coupling of this sudden change to 
numerical time derivatives drives instabilities in our code, 
in a manner reminiscent of the initial burst of radiation 
that occurs due to the sudden appearance of the particle 
at the start of our evolution; see Fig. 2 of Paper I. 

This problem is substantially mitigated by using a 
delta representation with a wider stencil; examples of 
this are described in Paper I. Wide stencils reduce the 
amount by which each weight changes from step to step, 
thereby reducing numerical noise. Another useful tool is 
to increase the order of the delta representation, thereby 
increasing the smoothness of the delta and its derivatives. 
This is particularly important since the Teukolsky equa- 
tion is a second-order differential equation; some smooth- 
ness in the derivatives is necessary to prevent the differen- 
tial operator from seeding excessive noise. Finally, resid- 
ual high frequency noise can be removed by convolving 



the source with a low pass filter 2 . These three techniques 
are each described in the following subsections. 

Each of these techniques smear out the delta function, 
pushing us away from the idealization of a zero width 
singularity. Choosing between stability (which tends to 
push us to a wider delta) and faithful representation of 
the singularity (which pushes us to a narrow delta) leads 
us to an optimization problem; we tune our delta repre- 
sentation in a way that (hopefully) minimizes numerical 
noise and maximizes accuracy. Note also that, in addi- 
tion to high-frequency noise generated by abrupt move- 
ment of the delta across the grid, spurious excitations of 
the quasinormal modes of the black hole also appear due 
this motion. This source of "noise" appears to be con- 
trolled by grid resolution — wider grids lead to less point- 
like deltas, which spuriously excite these modes. This 
spurious contribution to the EMRI waves can be miti- 
gated with a form of Richardson extrapolation [2Cj . We 
discuss this further in Sec. IHII and the Appendix. 



A. Higher order delta functions 

Discrete delta representations based on linear and cu- 
bic interpolation were derived in Paper I. We now extend 
this process to arbitrary polynomial order, equipping us 
with an entire family of discrete delta functions. 

As in Paper I, we start from the defining integral, 



a(t)+e 



a(i)-6 



dx f(x) S[x - a(t)} = f[a(t)} . (2.7) 



Let Xk+n-i < oi < Xk+n] the reason for our somewhat 
idiosyncratic choice of subscripts will become clear as we 
proceed. For clarity, we will not explicitly write out the 
time dependence of a; the reader should bear in mind 
that a = a(t) in all that follows. Rewriting Eq. (|2.7p as 
a sum over a finite step size, we have 

/ dx f(x) S(x - a) ~ h^J (xj)5i 

Ja-e i 

~ h^2f(xi)Si . (2.8) 



The function /(a) can be approximated by the La- 
grange interpolating polynomial, 



/(«) 



fc+2n-l 

E 

i— k 



11(a) 



(a - Xi)U'{xi) 



f(xt) 



(2.9) 



2 An obvious brute force workaround left off this list is to simply 
make the grid extremely fine and use tiny time steps. This does 
not address the root cause of instabilities seeded by particle mo- 
tion, though it is certainly something used in practice (to the 
extent that computational limits allow). 
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where 2n is the order of interpolation and 

fc+2n~l 2n— 1 

n(a) = [| (a - Xi) = JJ (a - x fc+l ) (2.10) 

z= jfe i= 

= n (xj-xk). (2.ii) 



II' fo) = 



dn 



da 



Inserting this in Eq. (12. 8|) leaves us with 

fc+2n-l 



E 

i—k 



±±—f(xi) = h^2f(x i )S i ;(2.12) 



(a - Xi)W(xi) 



comparing coefficients of f(xi) allows us to read off 6. 

11(a) 



h(a - Xi)U'(xi) ' 



2.13) 



We thus see that Si is non-zero for i G [k, k + 2n — 1]. 

The weights for derivatives of the delta function can 
be obtained similarly. Writing the identities 



dx f(x) S'(x -a) = -f'(a) 
dx f(x) 6"(x - a) = f"(a) 



(2.14) 
(2.15) 



as sums gives us 

hj2 /(a*) ^ - -j» 



-hY^f{xi)& 

i 

fe+2n-l 



E /(**) 4 = - E 



IL(a)f{ Xi 



h(a - sct)II'(a;t) 



(2.16) 



hJ2 fi?i)6" - /» 

= hj2 f"(?i)&i 



fc+2n-l 



/(^) = e 



i—k 



n(*)f"(.xi) 

h(a - x l )W(x i ) 



(2.17) 



We now insert centered finite difference formulae for the 
derivatives of f(xi) to obtain 



fc+2n-l 

E ^ = - E 



E(a) 



ft(a - a;i)II'(a:j) 



/(Xi+l) - / (Xj-l) 



2ft 



(2.18) 



fc+2n-l 

E /(**) - E 



H(a) 



ft(a - ari)n'(ii) 
f(x i+1 ) - 2f(xi) + f(xi-i) 



ft 2 



(2.19) 



Expressions (|2. 18|) and (|2. 19[1 are in a form that makes 
it simple to read off 8[ and 8". For example, Sj can be 
calculated by setting f(xj) = 1 and f(xi) — 0, I ^ j. 
It is straightforward to verify that setting n = 1 and 
ri = 2 reproduces the weights given by the two-point 
linear hat and the cubic formulae (described in Paper 
I) respectively. We also note that the delta derivative 
coefficients are non-zero for i £ [k — 1, k + 2n]. 



B. Wider stencils at a given interpolation order 

In Paper I, we generalized the two-point linear hat 
delta function such that it can be represented over a 
larger number of points. Similarly, we develop a pro- 
cedure to widen the stencil of the generalized model ob- 
tained from Eqs. (|2~T5j) . (|2~T5|) . and (JHHJ). 

Consider a model for 8i obtained from Eq. (|2.13[) for 
some n = m. Then, Si ^ for i E [k, . . . , k + 2m — 1]. 
Our goal is to widen this representation by some integer 
factor w such that the coefficients are non-zero for a wider 
range of grid points. Let us label the weights of this wider 
representation by Sf, with Sf ^ for i G [k,...,k + 
2wm — 1] . It should be emphasized that this is different 
from simply using Eq. (|2.13[) with n = wm; we have not 
changed the polynomial order, it remains fixed at 2m. 

For concreteness, let us choose w = 2, doubling the 
number of points in the delta representation. We infer 
the coefficients Sf at gridpoints i = k, k + 2, k + 4, . . ., 
k + Am — 2, by widening the grid by a factor of two: We 
evaluate Si with ft — > 2ft, Xk+j — ► Xk+2j to get 



$k+2j — 3k+j\h^2h,x k+j ->x k+2j 

11(a) 



where 



2ft(a - x k+2j )Tl'(xk+2 3 ) 



n(a) = (a - x k+2 i) 

i=0 



(2.20) 



(2.21) 



Finally, we need Sf at the intermediate points i = k + 1, 
fc + 3, . . ., fc + 4m — 1. We do this by exploiting the trans- 
lational symmetry of the problem. Momentarily reinsert 
the time dependence of the S's and a. Now consider the 
hypothetical situation where 



a (to) = «o , 
a(ti) = ao — ft 



(2.22) 



i.e, a(t) changes by a grid spacing from to to t\. We must 
have 



<*k+2j(*l) 



^ 2 



k+2j 



(to)Ja(t)- 



■a — h 



<>k+2j+l(~ t o) 

S k+2j+l (*0) J a(t)^cJ 2 - 23 ) 



We can turn this equation the other way around to read 
off the coefficient 8f +2 j+i a * to'- Simply replace a(t) with 
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a(t) — h in the formula for 5f. +2 j{to) to obtain S%. +2 j +1 (to). 
Since there was nothing special about our time slice, to, 
we find 

<*fe+2i+l(*n) = S k+2 3 ( t n)\a(t)^a(t)-h (2-24) 

for any moment t n . 

Though we chose w = 2 for concreteness, the above 
argument can be generalized to any integer w. Since our 
result holds for all time slices, we again suppress the time 
dependence to obtain expressions for any integer w: 



2m- 1 



11(a) = JY (a - Xk+wi) 



»=o 



r«) r I 

v k+wj ~ u k+jjh^wh,x k+j ^x k+u 

11(a) 



(2.25) 
(2.26) 
, (2.27) 



tuft(a - x k+wj )W(x k+wj ) 

ft" Jf2 I 

°k+wj+l — °k+wjJa(t)—>a(t)-lh 

for / G [1,2,. ..,to - 1] . (2.28) 

These techniques carry over to the derivatives as well: 

^k+wj = ^k+j\h—rwh,x k+ j—>-x k+wi (2.29) 

slw j:fw I 

°k+wj+l — °k+wj\a(t)^a(t)-lh 

for I G [l,2,...,w-l] ; (2.30) 



and 



°k+wj 



J fe+j'J h^wh,x k+J ^x k+mj (2-31) 
"fc+ujj+Z — °fc-HojJa(t)-«3((i)-I/i 

for/e [l,2,...,w-l] . (2.32) 



These should be used with Eqs. (|233]l . (|2~18|) and (|2~T9|) 
to widen the Teukolsky source term by any factor w. 



On a discrete grid, this becomes 

s(xi) = Si = gi(a)6i + g2(a)S' i + g^{a)5'l . (2.36) 

If the delta function and its derivatives span 2n + 2 grid 
points, with Xk+ n -i < a < %k+n, then s, ^ for i G 
[fc — 1, . . . , fe + 2n]. The source s, is zero everywhere else 
on the grid. 

The Gaussian filter is given by 



Ck 



exp[- {kh/bf /2] 
EL- P exph (*V^) 2 /2] 



(2.37) 



where fc G [— p, — p + 1, . . . ,p] and b is the width of the 
filter. The quantities p and b are adjustable parameters. 
Typically, we use p = 30 and 6 = 1.5/i. Notice that 



E 

i=—p 



Ci 



(2.38) 



this normalization guarantees that the integrated value 
of any function convolved with the filter is unchanged. 
We now convolve the source with the filter to obtain 



sgk 



Sk+i 



(2.39) 



where sgk is the smoothed source term. This indicates 
that sgk ^ for k G [fe — p, . . . , k + 2n +p — 1]. 

A wide filter spreads the source over a large domain on 
the numerical grid and thus increases errors, although it 
eliminates spurious harmonics. We have found that using 
a wide stencil followed by a narrow Gaussian smoother 
works very well to reduce numerical noise and minimize 
errors from an insufficiently pointlike source. 



C. Smoothing the source with a Gaussian filter 

Further control of numerical noise can be achieved by 
filtering high frequency components in the source term. 
This requires a convolution of the source with a discrete 
low pass filter. We use a Gaussian filter because it max- 
imizes the uncertainty principle — it can be localized in 
both position and frequency with greatest efficiency. 

Consider a source of the form 

s(x) = f 1 (x)5{x-a) + f 2 (x)5'(x-a) + f 3 (x)5"{x-a) . 

(2.33) 

Delta function identities allow us to rewrite this as 

s(x) — gi(a)S(x — a) + g2(oi)S'(x — a) + g^(a)5" (x — a) , 

(2.34) 

where 

9l (a) = h(a)-f'(a) + tt(a) , 

52(a) = f 2 (a)-2f 3 (a), 

93(a) = fS(a) . (2.35) 



D. Order of convergence of the filtered delta 

Paper I discussed in detail the convergence of a code 
that uses a discrete delta. Crucial background is given 
by Ref. [2l| and summarized in Paper I. The key point 
is that the moment 

fc+2n-l 

M r = h S^x.-aY (2.40) 

i—k 

controls the delta's convergence properties. Clearly, 
Mq = 1 (otherwise the delta is not properly normalized); 
in the continuum limit, M r = for r > 0. For the dis- 
crete delta, the smallest non-zero value of r for which 
M r ^ sets the order of convergence. In particular, if 
M r ^ 0, then a code which uses this delta will be no 
higher than rth-order convergent. 

We now show that, if a delta representation is second- 
order convergent before smoothing with the Gaussian fil- 
ter (Mo = 1, Mi = 0, M2 7^ 0), it will remain second- 
order convergent after smoothing. Upon convolving the 
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discrete delta with the Gaussian smoother, we find 

p 



3 = ~P 



(2.41) 



Let us denote the moments of the smoothed delta by M-!. 
As discussed in Sec. Ill CI the convolution does not change 
the delta's normalization as long as the Gaussian filter is 



itself properly normalized; thus 
fc+2n-l 

M$ = h 5 9i = 1 



(2.42) 



i—k 



We now examine the next higher moment of the 
smoothed delta: 



fc+2n+p-l 



Mf = /i 2^ Sgi(xi - a) = h ^ ^ c^i+j^ - a) , 



i=k—p j—^P i—k—p 

p fc+2n+p— 1 

j——p i—k—p 
p k-\-2n-\-p— 1 

= ft X C J X <5< +J - (a: i+ j - a - j/i) , 

j——p i—k—p 

p k-\-2n-\-p— 1 p fc+2n+p— 1 

J=— P i=k~p j=— p i=k—p 

The first term on the final line of (|2.43[) gives zero: Since &l x l = a > 

p k+2n-\-p— 1 p fe+2n+p+j— 1 

^ X C J X 4fifa+j -a) = ft X C J X Sjfa-a) 

j'=— P i=k-p j=-p l=k-p+j 

= . (2.44) 

The second line follows because \j\ <= p, 5i = if i lies outside [k, k + 2n — 1] and Mj = 0. 
The second term on the final line of (|2.43[) also yields zero: 

p fe+2n+p-l p k+2n+p+j-l 

ft X h i°j X Si +j = h2 X ^ X 5 ' ' 

j=-p i=k-p j=-p l=k-p+j 

P 

= ft 2 X ^ 

3=—P 

= . (2.45) 



The Gaussian filter's symmetry property Cj = C-j has been applied in the last step. Hence, we find Mf = Mi = 0. 



Evaluating the second moment proceeds similarly, but 
we find in the end terms involving *YTj = _ p j 2 Cj which do 
not vanish. Thus, M| is the first non- vanishing moment 
of the discrete delta, demonstrating that the Gaussian- 
filtered discrete delta function exhibits second-order con- 
vergence. The argument can be extended to the delta 
derivatives as well. The smoothed Teukolsky source term 
will thus be second-order convergent. 



III. WAVEFORMS AND COMPARISONS FOR 
GENERIC GEODESIC KERR ORBITS 



We now present the waveforms generated by a point 
particle in a geodesic orbit around a Kerr black hole. 
The code used to generate these waves is discussed in 
detail in Paper I; the only important change to that dis- 
cussion is that the source term uses the techniques pre- 
sented in Sec. above. We begin by reviewing Kerr black 
hole geodesies, sketching the numerical scheme used to 
solve the equations of motion. We then examine differ- 
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ent classes of eccentric and inclined orbits and compare 
the waveforms against those obtained from a frequency- 
domain code whose details are given in Ref. lj. We 
compute the correlation between the two waveforms in 
order to measure our level of agreement with frequency- 
domain waveforms. 

Our numerical grid is laid out in Boyer-Lindquist co- 
ordinats and uses (6r,66,5t) = (0.04M, tt/60, 0.02M) for 
the radial, angular and temporal resolutions. The source 
term is constructed using Eqs. ([235]) . (gHHD and (HS) 
with ng in the range 3-9 (depending on the orbit) for the 
angular delta-function and n r = 2 for the radial delta. 
We use a Gaussian filter of width b = 1.559 to smooth 
higher harmonic noise. 



eliminate the turning points by remapping the coordi- 
nates r and to parameters which accumulate secularly. 
The following parameterization, inspired by the Newto- 
nian limit, has been found to work extremely well even 
deep in the strong field of rapidly rotating black holes: 



P 



1 + e cos tjj 



cos 9 — cos# min cosx 



(3.5) 
(3.6) 



In the Newtonian limit, p is the orbit's semi-latus rectum, 
and e is its eccentricity; 9 m i n is the minimum value of 9 
reached by the orbiting body, and is used to define the 
orbit's inclination 9;„ r 



A. Geodesies in Kerr spacetime 

The source term for the time-domain code takes as 
input the worldlinc of the perturbation's source. Here, 
we neglect radiation reaction and assume that the point 
particle follows a bound geodesic trajectory around the 
central massive black hole. This bound trajectory can 
be computed by numerically integrating the geodesic 
equations. We now briefly review how we massage the 
geodesic equations to put them into a form that makes 
for accurate numerical calculation; this material is pre- 
sented in greater depth in Sec. IIC of Ref. 

The normal "textbook" presentation of the equations 
governing Kerr black hole geodesies is 



drY 



d9Y 



,dt 



[E(r 2 +a 2 )- 

-A [r 2 + (L z 
R(r) 



aL z ] 
-aE) 2 



Q] 



(3-1) 



cos 2 9 [a 2 (l - E 2 ) + L 2 J sin 2 9] 



(3.2) 



sin 



~aE+— [E(r 2 + a 2 ) - aL z ] 



a{L z - aE sin 2 9) 



[E(r 2 + a 2 ) - aL z ] 



[See, e.g., Ref. (H, Eqs. (33.32a-d).] Here, £ 



(3.3) 



(3.4) 



a 2 cos 2 9, A 



2Mr + a 2 (where a = \S\/M is the 



black hole's spin angular momentum per unit mass). The 
constants of motion are orbital energy E, axial angular 
momentum L z , and Carter constant Q. 

This form of the equations of motion is not well suited 
to numerical studies; in particular, dr/dr and d9 jdr pass 
through zero and change sign when the orbiting body 
goes through turning points associated with those mo- 
tions. A handy way to eliminate these problems is to 



(3.7) 



Once E, L z , and Q are specified, p, e, and 9[ nc are 
fully determined. It is then a straightforward matter to 
turn Eqs. (|3.ip and (|3.2[) into expressions for dip/dr and 
dx/dr; see Ref. [19( for details. The resulting expres- 
sions behave extremely well for all bound orbits outside 
the black hole's event horizon. A numerical integrator 
for these variables allows us to compute the dynamics of 
our orbiting body's Teukolsky equation source term. 

Before moving on, we note that, within the context of 
the dissipative-only or radiative approximation to inspi- 
ral, it is simple to modify these equations to build the 
worldline of an inspir ailing body: We simply allow the 
orbital "constants" (E, L z , and Q; or, p, e, and 6*j nc ) to 
evolve according to the inspiral law. Reference [lj| uses 
approximate radiation reaction, based on fits to strong- 
field radiation reaction calculations in regimes where it is 
well understood, to compute the inspiral worldlines which 
underlie the "kludge" waveforms. We use this prescrip- 
tion for evolving the constants in Sec. lIVI to demonstrate 
this code's ability to compute inspiral waves. 



B. Comparison with frequency-domain waveforms 

To validate our waveforms, we compare with the "snap- 
shots" generated using the frequency-domain code de- 
scribed in Ref. This code uses the fact that bound 
Kerr geodesies are fully described by three frequencies 
(radial fi r , latitudinal Qg, and axial f^) to build the 
waveform from a geodesic orbit as a sum over harmonics 
of these frequencies [231 ] . Since both the time-domain and 
frequency-domain codes solve the same master equation, 
they should produce identical waveforms for identical or- 
bits, so long as each code is sufficiently accurate. 

To quantify the accuracy with which a time-domain 
waveform X = (xi,X2, ■ ■ ■ , x n ) agrees with a frequency- 
domain waveform Y — (yi, y^, ■ ■ ■ , y n )i we use the follow- 
ing correlation measure: 



rxY 



- x)(vi - v) 



wuvi-vy 



(3.8) 
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(The sums in all cases are from i = 1 to i = n.) This co- 
efficient is identical to the match between two waveforms 
defined by Owen [24| in the white noise limit [noise spec- 
tral density Sh(f) — constant]. One might expect the 
waveforms' mean values x and y to equal zero. However, 
finite duration effects can make these quantities slightly 
non-zero, so it is useful to explicitly do this subtraction. 

A useful reformulation of Eq. (|3.8p is 



rxv 



> E x iVi "S^E!/! 



(3.9) 



Note that txy is always between —1 and 1; a value close 
to 1 indicates that the two waveforms are well corre- 
lated. Note also that the correlation depends on how 
many points n are used in comparing the two waveforms 
(or equivalently, the span of time over which we compare 
the waves). We have found that as long as n > sev- 
eral hundred, we get consistent results: Changing n for 
a given comparison only causes small variations in the 
fourth significant digit of rxY ■ 

It is of course possible to concoct other measures of 
how well two waveforms agree. Ideally, disagreements 
between waveforms should be quantified in terms of their 
observational significance. For example, Cutler and Val- 
lisneri have demonstrated that it is not unusual for wave- 
forms with a match of 0.9999 to differ significantly in 
their estimates of the parameters which describe the 
source 25]. For our present purpose, txy is sufficient to 
demonstrate that our time-domain code produces high 
quality waveforms; whether they are sufficiently high 
quality to be used for GW measurement purposes will 
need to be re-examined at a later time. 

An important step in producing accurate waveforms 
is to perform runs at multiple resolutions, then esti- 
mate (and eliminate) the waveform error using a form 
of Richardson extrapolation [2(|. This plays a crucial 
role in reducing "noise" from spurious excitation of the 
large black hole's quasinormal modes. The details of this 
extrapolation technique are described in Appendix [A"l 

Tables ffl El EH El El and ED list the correlation co- 
efficients for the m — 2 and m = 3 azimuthal modes of 
different classes of orbits. The coefficient is greater than 
0.99 for a large fraction of parameter space. Time do- 
main runs corresponding to each column required about 
125 CPU hours on an Apple MacPro processor. That 
code was compiled using the Intel CH — h compiler. The 
frequency domain code's cost is about 3-4 CPU hours 
per waveform when attempting to get both asymptotic 
energy fluxes to accuracies of about 0.1% to 1% on a ma- 
chine using a 3.2 GHz Pentium 4 Xeon processor. We 
also show (Figs. [U [21 and[3|) examples of the waves, com- 
puted with both time- and frequency-domain codes, to 
give the reader a visual sense of the overlap. 



TABLE I: Correlation between time- and frequency-domain 
waveforms for the m — 2 mode for a range of equatorial, ec- 
centric orbits. The parameters p, e and #j nc are semi-latus rec- 
tum, eccentricity, and inclination of the geodesic orbit, a/M 
is the black hole spin and 6d is the angle between the spin axis 
and the line of sight to the observer. The last two columns 
show correlations for the plus and cross polarizations. 



p/M 


e 


ft- l At±cr\ 
fine (, Cleg J 


a i ivi 


Vd (aeg) 


corr. 


h x corr . 


0.4 ( z 


n q 
U.o 


U 


U.o 


QO 


n not; i 

u.yyoi 


11.9902 


6.472 


0.3 





0.3 


60 


0.9969 


0.9969 


6.472 


0.3 





0.3 


90 


0.9974 


0.9975 


5.768 


0.3 





0.7 


30 


0.9971 


0.9971 


5.768 


0.3 





0.7 


60 


0.9977 


0.9978 


5.768 


0.3 





0.7 


90 


0.9983 


0.9983 


6.472 


0.7 





0.3 


30 


0.9915 


0.9911 


6.472 


0.7 





0.3 


60 


0.9911 


0.9908 


6.472 


0.7 





0.3 


90 


0.9900 


0.9901 


5.768 


0.7 





0.7 


30 


0.9625 


0.9607 


5.768 


0.7 





0.7 


60 


0.9621 


0.9601 


5.768 


0.7 





0.7 


90 


0.9596 


0.9578 



TABLE II: Correlation between time- and frequency-domain 
waveforms for the m — 2 mode for a range of inclined nearly 
circular orbits. All symbols have the same meaning as in 
Table U 



p/M 


e 


fine (deg) 


a/M 


Od (deg) 


h+ corr. 


h x corr. 


6 


10" 4 


45 


0.5 


60 


0.9968 


0.9967 


6 


10" 4 


45 


0.5 


90 


0.9961 


0.9960 


8 


10" 4 


45 


0.5 


60 


0.9923 


0.9919 


8 


10" 4 


45 


0.5 


90 


0.9908 


0.9903 


6 


10" 4 


45 


0.9 


60 


0.9967 


0.9967 


6 


10" 4 


45 


0.9 


90 


0.9961 


0.9961 


8 


10" 4 


45 


0.9 


60 


0.9920 


0.9919 


8 


10" 4 


45 


0.9 


90 


0.9905 


0.9907 


6 


10" 4 


60 


0.5 


60 


0.9964 


0.9965 


6 


10" 4 


60 


0.5 


90 


0.9952 


0.9952 


8 


10" 4 


60 


0.5 


60 


0.9917 


0.9910 


8 


10" 4 


60 


0.5 


90 


0.9888 


0.9882 


6 


10" 4 


60 


0.9 


60 


0.9986 


0.9986 


6 


KT 4 


60 


0.9 


90 


0.9981 


0.9982 


8 


10" 4 


60 


0.9 


60 


0.9917 


0.9915 


8 


10" 4 


60 


0.9 


90 


0.9891 


0.9890 



IV. INSPIRAL WAVEFORMS 

Having demonstrated that the finite-impulse source 
works well for astrophysically relevant generic black hole 
orbits, we now examine how well we do evolving through 
a sequence of such orbits. Since each orbit in the sequence 
is no different than the orbits that we validated against 
in Sec. IIIIB1 we anticipate no great difficulty here. In- 
deed, the biggest challenge is choosing a method to evolve 
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Time and frequency domain waveforms extracted at 8 = jt/2. 
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time/M 
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time/M 
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time/M 

FIG. 1: Comparison of time- and frequency-domain waveforms. We show waves for the m = 2 mode from a point particle with 
orbital parameters p — 6.472M, e = 0.3 and Sine = orbiting a black hole with spin a/M = 0.3. The angle between the spin 
axis of the black hole and the line of sight is 8d — tt/2. Time-domain results are in black, frequency-domain results in red. Top 
panel: "plus" polarizations in dimensionless units. Middle: "cross" polarizations. Bottom: Comparison of \h+ — ihx\- This 
last quantity gives a good visual measure of the level of agreement between the two waveforms. The correlations between the 
two waveforms are 0.9974 (plus) and 0.9975 (cross). 



TABLE III: Correlation between time- and frequency-domain 
waveforms for the m — 2 mode for a range of generic orbits. 
All symbols have the same meaning as in Table [T] 



p/M 
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Sine (deg) 


a/M 


0d (deg) 


h+ corr. 


h x corr. 
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0.3 


40 


0.9 


60 


0.9978 


0.9978 
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0.3 


40 


0.9 


90 


0.9976 


0.9976 
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40 
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60 


0.9898 


0.9897 
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0.3 


40 


0.5 


90 


0.9910 


0.9910 


6 


0.7 


40 


0.9 


60 


0.9898 


0.9906 


6 


0.7 


40 


0.9 


90 


0.9889 


0.9891 


6 


0.7 


60 


0.9 


60 


0.9905 


0.9868 


6 


0.7 


60 


0.9 


90 


0.9895 


0.9866 


6 


0.3 


60 


0.9 


60 


0.9961 


0.9962 


6 


0.3 


60 


0.9 


90 


0.9950 


0.9954 


8 


0.3 


60 


0.5 


60 


0.9906 


0.9890 


8 


0.3 


60 


0.5 


90 


0.9884 


0.9866 



TABLE IV: Correlation between time- and frequency-domain 
waveforms for the m = 3 mode for a range of equatorial ec- 
centric orbits. All symbols are as in Table [I] 



p/M 


e 


6> inc (deg) 


a/M 


0d (deg) 


h+ corr. 


h x corr. 


6.472 


0.3 





0.3 


30 


0.9908 


0.9909 


6.472 


0.3 





0.3 


60 


0.9922 


0.9922 


6.472 


0.3 





0.3 


90 


0.9930 


0.9931 


5.768 


0.3 





0.7 


30 


0.9934 


0.9935 


5.768 


0.3 





0.7 


60 


0.9943 


0.9944 


5.768 


0.3 





0.7 


90 


0.9948 


0.9948 


6.472 


0.7 





0.3 


30 


0.9931 


0.9931 


6.472 


0.7 





0.3 


60 


0.9905 


0.9906 


6.472 


0.7 





0.3 


90 


0.9923 


0.9923 


5.768 


0.7 





0.7 


30 


0.9928 


0.9929 


5.768 


0.7 





0.7 


60 


0.9932 


0.9930 


5.768 


0.7 





0.7 


90 


0.9920 


0.9921 



through our sequence. Our goal is to do this with a 
frequency-domain code to build the orbital-constant tra- 
jectory [E(t), L z (t),Q(t)]. To quickly produce results 
that are qualitatively correct, we presently make this tra- 



jectory using the "kludge" inspiral treatment described 
in Ref. [HI], and used to make model waveforms in Ref. 
[l9l |. The "kludge" uses a somewhat idiosyncratic mix 
of post-Newtonian backreaction formulae combined with 
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FIG. 2: Comparison of time- and frequency-domain waveforms. Here, we show waves for the m = 2 mode for a geodesic with 
p = 6A/, e = 0.3 and 6>i n c = 7r/3 about a black hole with spin a/M = 0.9; black is time-domain results, red is frequency domain. 
The correlations in this case are 0.9961 (plus) and 0.9962 (cross). 



numerical results from frequency-domain backreaction in 
the circular, inclined (e = 0, f?i nc ^ 0) and eccentric, 
equatorial (e ^ 0, 9 mc = 0) limits to estimate the prop- 
erties of EMRI waves. By construction, the results agree 
very well with Teukolsky-based inspirals in those limits; 
for the generic case, they produce plausible inspirals. 

Figure [4] shows our waveform for a "kludge" inspiral. 
We took the large black hole to have spin a = 0.5M, and 
set the mass ratio to n/M = 0.016. The orbit was ini- 
tially chosen to have semi-latus rectum p = 10M, eccen- 
tricity e = 0.5, and inclination f?i nc = 0.5 radians. This 
figure shows features reminiscent of the geodesic snap- 
shots shown in Figs. [TJ [H and [3J in addition, one can 
clearly see evolution of the wave's properties. The in- 
crease in the wave's frequency, largely due to the decay of 
the orbit's semi-latus rectum, is quite clear. Perhaps less 
obvious is a signature of the eccentricity's decay. This 
is illustrated most clearly by comparing the lower left 
and lower right panels of Fig. [H which zoom onto early 
and late portions of the inspiral. Early on, the wave- 
form is dominated by a series of high-frequency bursts; 
these occur when the small body passes through periapsis 
and "whirls" most rapidly about the massive black hole. 
There is then a relatively quiet section while the body 
"zooms" out to apoapsis, and then comes in to "whirl" 
at periapsis again. As eccentricity shrinks, the difference 
between periapsis and apoapsis becomes smaller. The 



high-frequency bursts crowd closer and closer together, 
approaching a continuum sinusoid as the eccentricity ap- 
proaches zero. 

Although this inspiral model is somewhat unphysical, 
we expect that it shares many properties with true adia- 
batic inspiral waveforms. In particular, the spectral evo- 
lution of a wave like that in Fig. 2] should be quite sim- 
ilar to the evolution of real EMRI waveforms. It should 
be emphasized that computing the waveform shown in 
Fig. 0] required about as much computational effort as 
computing the geodesic snapshot waves, Figs. [U [21 and 
[3] (modulo a factor ~ 4-5 since the waveform in Fig. Q] 
lasts ~ 4-5 times longer than the others) . Given a robust 
code to generate the inspiral worldlinc of EMRI systems, 
the waveforms that our code produces should be a useful 
tool for examining issues in LISA measurement and data 
analysis. 



V. SUMMARY AND FUTURE WORK 

We have now shown that the finite impulse delta repre- 
sentation of the time-domain Teukolsky equation's source 
works very well for complicated and astrophysically rel- 
evant orbits. In our previous analysis (|ll|. Paper I), 
we confined ourselves to the simplest circular, equato- 
rial black hole orbits. The basic ideas from Paper I 
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FIG. 3: Comparison of time- and frequency-domain waveforms. These waves are for the m = 3 mode from a circular geodesic 
with orbital parameters p = 6M, and 8i nc — 7r/4 around a hole with spin a/M = 0.9. All symbols have the same meaning as 
in Fig. Q] The correlations are 0.9769 (plus) and 0.9770 (cross). 



work well even when the source arises from highly in- 
clined and highly eccentric orbits, and when the source 
evolves through a sequence of those orbits. It is now a 
relatively straightforward matter to compute the waves 
arising from a body following any reasonably behaved 
worldline in the spacetime of a black hole. 

The primary complication arising from these more 
generic orbit classes is that the orbiting body will cross 
zones within our numerical grid. The source thus be- 
comes dynamical; the finite-impulse delta must likewise 
be dynamical to represent it. The evolution of the im- 
pulses that we use to represent the delta can seed numer- 
ical noise, reducing the calculation's accuracy. We have 
found that minor extensions of Paper I's basic techniques 
greatly mitigate the impact of this source of numerical 
noise. In particular, by using a higher-order representa- 
tion (Sec. Ill Ap . the delta is smoothed enough that the 
coupling to the Teukolsky equation's second-order differ- 
ential operators does not seed much error. Widening the 
delta's stencil (Sec. IIIBI) also helps, since the fractional 
change in a given impulse will be less if the delta is repre- 
sented by more impulses. Finally, residual high frequency 
noise not removed by these techniques can be taken out 
by convolving the Teukolsky source term with a low-pass 
(Gaussian) filter (Sec. Ill C|) . It's worth emphasizing that 
we smooth the entire source term, not just the delta func- 
tion (which would arguably make our delta rather similar 



to the truncated Gaussian [26|, l27| which this technique 
was designed to improve upon). 

Comparison with results from the frequency-domain 
[lij demonstrates that the waveforms generated with this 
source term are of very high quality (Sec. IIII[) . Visually, 
the waveforms lie on top of one another in every case 
that we have examined; a quantitative overlap integral 
demonstrates that waveforms from the two calculations 
are often more than 99% correlated. A key step in achiev- 
ing such high quality results is to estimate the largest er- 
rors in our time-domain calculations, and then subtract 
that estimate from our result. We do this by performing 
these calculations at two different grid resolutions; under 
the assumption that our dominant error is quadratic in 
grid spacing, we then estimate the magnitude of our er- 
ror (Appendix [X]). The excellent agreement we achieve 
with frequency-domain results validates this approach, at 
least for all the cases we have considered. 

So far, our main physics accomplishment is excel- 
lent agreement between time- and frequency-domain ap- 
proaches to waveform calculation. It should be empha- 
sized, however, that for waveform calculations, there will 
be a large set of circumstances in which time-domain 
codes are more efficient. For generic orbits, a frequency- 
domain code may require the calculation and summa- 
tion of many thousand multipoles and Fourier modes. A 
time-domain code "automatically" sums over all modes 
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FIG. 4: Waveform (m = 2 mode) of a small body spiraling into a massive black hole. We use "kludge" backreaction to evolve 
through a sequence of orbits, but compute the waves with our time-domain solver. The large black hole has spin a = 0.5M; 
the small body's orbit initially has parameters p — 10M, e = 0.5, and 9i nc = 0.5 radians. The mass ratio of the system is 
H/M = 0.016. The top panel shows the full span that we simulated; the bottom two panels are zooms on early (bottom left) 
and late (bottom right) segments. Note the clear evolution of the wave's frequency as the orbit's mean radius shrinks. 



(except the m index), so that (in principle) it is no more 
difficult to compute the waves from a highly inclined, 
highly eccentric black hole orbit than from an orbit with 
modest inclination and eccentricity. 

The real payoff of this tool will come when we allow the 
source to radiatively decay, evolving through a sequence 
of orbits. As a demonstration that this can be done, we 
use a "kludged" inspiral to compute a body's inspiral, 
and then use that inspiral as the source for our time- 
domain solver in Sec. IIVI Though not a physically accu- 
rate inspiral, this scenario shares many properties with 
the actual adiabatic inspiral. In particular, it demon- 
strates the computational advantage of a robust time- 
domain code for computing inspiral waveforms, given the 
worldline the inspiraling body follows. 

Future work will address our goal of complete wave- 
forms for the EMRI problem, in the context of the 
dissipation-only approximation to EMRI dynamics. We 
have recently extended our frequency-domain code to in- 
clude the evolution of Carter's constant in the radiative 



backreaction limit 15], and will use this code to pro- 
duce the radiation reaction data describing an inspiraling 
body. With this step in hand, no issue of principle stands 
in the way of coupling the time- and frequency-domain 
approaches to make usefully accurate EMRI waveforms. 
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TABLE V: Correlation between time- and frequency-domain 
waveforms for the m — 3 mode for a range of inclined nearly 
circular orbits. All symbols are as in Table [I] 
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#inc (deg) 
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0d (deg) 
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45 
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6 


itr 4 


60 


0.9 


60 


0.9665 


0.9661 


6 


itr 4 


60 


0.9 


90 


0.9680 


0.9678 


8 


itr 4 


60 


0.9 


60 


0.9463 


0.9473 


8 


itr 4 


60 


0.9 


90 


0.9608 


0.9641 



TABLE VI: Correlation between time- and frequency-domain 
waveforms for the m — 3 mode for a range of generic orbits. 
All symbols are as in Table [I] 



p/M 


e 


ft™ (deg) 


a/M 


0d (deg) 


h+ corr. 


hx corr. 


6 


0.3 


40 


0.9 


60 


0.9917 


0.9916 


6 


0.3 


40 


0.9 


90 


0.9915 


0.9914 


8 


0.3 


40 


0.5 


60 


0.9801 


0.9803 


8 


0.3 


40 


0.5 


90 


0.9785 


0.9785 


6 


0.7 


40 


0.9 


60 


0.9906 


0.9981 


6 


0.7 


40 


0.9 


90 


0.9899 


0.9895 


6 


0.7 


60 


0.9 


60 


0.9862 


0.9862 


6 


0.7 


60 


0.9 


90 


0.9819 


0.9821 


6 


0.3 


60 


0.9 


60 


0.9790 


0.9788 


6 


0.3 


60 


0.9 


90 


0.9840 


0.9839 


8 


0.3 


60 


0.5 


60 


0.9788 


0.9791 


8 


0.3 


60 


0.5 


90 


0.9747 


0.9744 



confirm the production results presented here. S. D.'s 
contribution to this analysis was carried out at the Jet 
Propulsion Laboratory, California Institute of Technol- 
ogy, under a contract with the National Aeronautics and 
Space Administration and funded through the internal 
Human Resources Development Fund Initiative and the 
LISA Mission Science Office. Some of the supercomput- 
ers used in this analysis were provided by funding from 
the JPL Office of the Chief Information Officer. 



APPENDIX A: WAVEFORM EXTRAPOLATION 

Here we describe the variation of Richardson extrapo- 
lation which we use to estimate and eliminate the largest 
errors arising from our finite difference scheme. In Ref. 
(llj . we showed that our algorithm is second order con- 
vergent. This means that we can write the solution at 
any given resolution as 



* t + a x 6r 2 + a 2 60 2 + a 3 SrS0 + 0(6 3 ) , (Al) 



where \& c is the computed solution and *f? t is the "true" 
solution. The final term 0(S 3 ) indicates that additional 
error terms will be third order in the grid spacing (and 
higher). The spatial and temporal dependences of ^f c 
and have been suppressed. We now perform runs at 
two different resolutions, (Sri,69i) and (5r2,802), with 
all other parameters fixed. The resolutions are chosen 
such that 



Sri 
Sr 2 



6di 

w 2 



n . 



(A2) 



Neglecting higher order terms, the two results can be 
written 



(A3) 
(A4) 

!, al- 



* c i ~ * t + ai5r\ + a 2 60l + a 3 6ri66i 
* c2 ~ * t + a x 6r\ + a 2 50 2 + a 3 Sr 2 60 2 

The relation between the two resolutions, Eq. 
lows us to write 



* c2 = *t + l/n 2 ( ai Sr 2 + a 2 S9 2 + a 3 Sn60i) . (A5) 
Subtracting Eq. (|A5|) from Eq. (|A3[) leaves us with 
*d - * c2 = (1 - l/n 2 )( ai 6r 2 + a 2 89\ + a 5 6n66i) ; 



rearranging, we find 



(axSrl + a 2 59\ + a^nSOi) = 



gel - ^c2 

1 - 1/n 2 



(A6) 



(A7) 



To the extent that neglect of higher-order errors is war- 
ranted, this estimates the largest source of error. Using 
Eq. (|A3|) we can now estimate the "true" value: 



(aiSrj + a 2 60\ 

gel ~ frc2 

1-1/n 2 ' 



a 3 6ri60i 



(A8) 



Figure [5] illustrates the improvement that this variant 
of Richardson extrapolation can yield. We plot h+ at 
two different resolutions: (6ri,66i) = (0.04, 7r/60) and 
(6r 2 ,80 2 ) = (0.026667, tt/90). We also show the ex- 
trapolated waveform, and the frequency-domain predic- 
tion. The particle is in a geodesic orbit with parameters 
p = 6M, 6 inc = 45°, e = 10~ 4 and the black hole has 
a spin of a — 0.9M. The two time-domain calculations 
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FIG. 5: Extrapolation applied to h + for the m = 3 mode 
from a point particle in a nearly circular geodesic with or- 
bital parameters e — 10~ 4 , p = 6M, and 8i nc = 7r/4 around 
a rotating black hole with spin a/M — 0.9. The dashed 
and solid black lines denote h+ obtained with resolutions 
(6r,S0) = (0.04,?r/60) and (0.026667, tt/90) respectively. The 
solid red line is the extrapolated waveform; the solid green 
line is the equivalent frequency-domain waveform. Notice 
how well the extrapolated time-domain wave agrees with the 
frequency-domain result (which is nearly hidden by the red 
curve) . 



each differ noticeably from the frequency-domain result; 
the extrapolated waveform by contrast agrees very well. 
This excellent agreement can be regarded as a modified 
three-level convergence test, whose first two levels are the 
time domain waveforms and third level is the frequency 
domain waveform. If the code were not second order con- 
vergent, our assumption for the functional form of the 
errors in Eq. (|A7[) would be erroneous. This would lead 
to a substantial disagreement between the extrapolated 
and frequency domain waveforms. 
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